Antiferromagnetic–ferromagnetic transition in zigzag graphene nanoribbons induced by substitutional doping
Yang Shenyuan1, 2, †, Li Jing1, Li Shu-Shen1, 3, 4
State Key Laboratory of Superlattices and Microstructures, Institute of Semiconductors, Chinese Academy of Sciences, Beijing 100083, China
School of Microelectronics, University of Chinese Academy of Sciences, Beijing 101408, China
College of Materials Science and Opto-Electronic Technology, University of Chinese Academy of Sciences, Beijing 101408, China
Synergetic Innovation Center of Quantum Information and Quantum Physics, University of Science and Technology of China, Hefei 230026, China

 

† Corresponding author. E-mail: syyang@semi.ac.cn

Project supported by the National Natural Science Foundation of China (Grant Nos. 11474274 and 61427901) and the National Basic Research Program of China (Grant No. 2014CB643902).

Abstract

Using first-principles calculations based on density functional theory, we show that the ground state of zigzag-edged graphene nanoribbons (ZGNRs) can be transformed from antiferromagnetic (AFM) order to ferromagnetic (FM) order by changing the substitutional sites of N or B dopants. This AFM–FM transition induced by substitutional sites is found to be a consequence of the competition between the edge and bulk states. The energy sequence of the edge and bulk states near the Fermi level is reversed in the AFM and FM configurations. When the dopant is substituted near the edge of the ribbon, the extra charge from the dopant is energetically favorable to occupy the edge states in AFM configuration. When the dopant is substituted near the center, the extra charge is energetically favorable to occupy the bulk states in FM configuration. Proper substrate with weak interaction is necessary to maintain the magnetic properties of the doped ZGNRs. Our study can serve as a guide to synthesize graphene nanostructures with stable FM order for future applications to spintronic devices.

1. Introduction

Magnetic carbon-based materials have attracted considerable research attention during the past decades.[1] Compared to the conventional d or f electron based magnetic systems, the magnetism of carbon-based materials is related to sp electrons, with the advantages of smaller spin–orbit coupling and longer spin coherent length. One of the interesting magnetic carbon-based materials is the graphene nanostructures with zigzag edges.[27] The carbon atoms on zigzag edge belong to the same sublattice of graphene, and will introduce flat bands near the Fermi level.[8] These flat bands correspond to localized edge states that are half-filled and spin-polarized, with parallel spins on the same edge.[2,4] In triangular zigzag-edged graphene nanoflakes, the spins on the three edges are ferromagnetically coupled due to topological frustration of the π-bonds, which is a generalization of the counting rule in magnetic organic molecules.[2] The magnetic coupling between two triangular nanoflakes that are linked by carbon chain can be ferromagnetic (FM) or antiferromagnetic (AFM), depending on the number of carbon atoms in the chain.[3] Compared to the triangular nanoflakes, one-dimensional zigzag-edged graphene nanoribbons (ZGNRs) are easier to synthesize via different methods.[57] In the ground state, the spins on the two edges of ZGNRs are predicted to be antiferromagnetically coupled, opening a bandgap in the otherwise metallic ribbons.[4,8] The opening of the bandgap and the edge states associated with the AFM coupling in ZGNRs have been confirmed by experiments, and the magnetic order is found to be stable up to room temperature.[6,7,9] Therefore, ZGNRs are considered as promising materials for future applications in graphene-based spintronic devices operating under ambient conditions.[6]

However, the AFM order in ZGNR results in zero magnetic moment and non-spin polarized transport, hindering its application to spintronic devices. Since the magnetism in ZGNRs is directly related to the spin-polarized edge states, various methods have been proposed to introduce non-zero magnetic moment by breaking the spin degeneracy of the edge states, such as substitutional doping,[1016] edge functionalization with chemical groups,[1719] geometric modification of the edges,[20] and external field.[21] Substitutional doping with heteroatoms such as N and B dopants is predicted to be stable at the edge,[22] which can obviously remove the spin degeneracy. Most previous studies have been focused on the edge substitution, and the spin-polarized transport has been predicted as expected.[1116] On the contrary, much less effort has been devoted to the far-edge substitution.[10,12,23] In most studies, an AFM coupling with broken inversion symmetry is usually supposed to be the ground state for the substitutionally doped ZGNRs, even though the energy difference between the AFM and FM states has not been carefully compared.[11]

ZGNRs with FM order are more desirable for spintronic devices. However, the FM order is less stable than the AFM order in ZGNRs. Several groups proposed different methods to convert the AFM ZGNRs into FM materials. Kan et al. found that the ground state of ZGNRs with line defects at the center can be transformed from AFM to FM by applying uniaxial strain or substituting with foreign atoms.[24] The transformation is ascribed to the defect mediated magnetic coupling, which can be tuned by strain or substitution. Nam et al. studied the magnetism of ground state of edge-terminated ZGNRs with organic radicals, and revealed a switch between AFM and FM configurations by additional N or B doping.[19] The switching of the magnetic order is explained by the spin alternation rule where adjacent carbon atoms show opposite spin directions. Another two groups reported a transition from AFM to FM configuration within a certain carrier doping range in ZGNRs, since the FM configuration do not have a bandgap and is more favorable by carrier doping.[25,26] However, the effects of substitutional sites on the magnetic properties in ZGNRs have not been carefully investigated.

In this paper, using first-principles calculations based on density functional theory, we show that the ground state of ZGNRs can be transformed from the AFM order to FM order by changing the substitutional doping sites. If N or B atom is doped near the edge of the ribbon, the ground state of the doped ZGNRs remains the AFM order, the same as the pristine ZGNRs. If N or B atom is doped near the center of the ribbon, the FM configuration of the doped ZGNRs becomes more stable than the AFM configuration. This AFM–FM transition induced by substitutional sites has never been reported before, and is found to be a consequence of the competition between the edge and bulk states. When the dopant is substituted near the edge, the extra charge introduced by the dopant is energetically more favorable to occupy the edge state in the AFM configuration. When the dopant is substituted near the center, the extra charge is energetically more favorable to occupy the bulk state in the FM configuration. We also discuss the effects of doping concentration and substrate on the magnetic properties of the doped ZGNRs. Our study might provide a guide to the experiment to synthesize graphene nanostructures with stable FM order, which is highly desirable for the future application of graphene-based spintronic devices.

2. Computational methods

We perform spin-polarized first-principles calculations using the Vienna ab-initio simulation package (VASP)[27] with the Ceperley–Alder local density approximation (LDA) for the exchange–correlation functional.[28] We employ the all-electron-like projector-augmented wave (PAW) pseudopotentials[29] and an energy cutoff of 450 eV for the plane-wave expansion. The Brillouin zone is sampled by Monkhorst–Pack mesh[30] of the form of 1 × 1 × K, with K = 6–30 depending on the length of ZGNR in the supercell along the growth direction. All the atomic coordinates are fully relaxed until the Hellmann–Feynman forces on atom are less than 0.001 eV/Å. Different initial magnetic couplings on the two edges are considered to get the energy difference between AFM and FM configurations for each doped ZGNR. We calculate selective systems using more accurate HSE06 functional, and find that LDA predicts the same structure stability between the AFM and FM configurations for the doped ZGNRs.

A one-dimensional periodic boundary condition is used along the growth direction (z direction) of ZGNRs (Fig. 1). A vacuum larger than 15 Å is introduced along both the x and y directions to avoid interactions between images. The edges of the ribbons are passivated by H atoms. We name the ZGNR with a ZM × N notation similar to Ref. [10], where M and N are the numbers of C lines across the ZGNR width and length in the supercell, respectively. The substitutional site for N and B dopants is denoted by numbers 1, 2, 3, … from the edge to the center, as shown in Fig. 1. The N and B doped ribbons are denoted by N-ZM × N and B-ZM × N, respectively.

Fig. 1. Scheme of one N atom substituting at the edge of Z6 × 6 zigzag nanoribbon. The substitutional sites are denoted by numbers 1–6, from the edge to the center of the nanoribbon. The gray, blue, and white balls represent the C, N, and H atoms, respectively.
3. Results and discussion

As predicted by previous calculations, the ZGNR has an AFM ground state with opposite spins on the two edges.[4,31,32] Figure 2(a) shows the band structure of pristine Z6 × 6. The charge density of the two flat bands near the Fermi level (EF) is mostly localized at the edges, and decays into the bulk of the ribbon. These states are referred to edge states hereafter. The spin-up and spin-down components of the flat bands are degenerate in energy, but separated in space, with spin-up states localized on one edge and spin-down states on the other edge. The AFM order of ZGNRs is consistent with Lieb’s theorem for bipartite system.[33] The atoms on the same edge belong to the same sublattice of graphene and have parallel spins, while the atoms on different edges belong to different sublattices and have opposite spins. The spin splitting of the edge states opens a band gap near k = 2π/3a (a is the lattice constant of graphene).[4,31] On the other hand, the charge density of the other dispersed bands is distributed across the whole ribbon, and thus is referred to bulk states. As seen from Fig. 2(a), in the AFM ground state, the occupied edge states (flat band below EF) have higher energies than the occupied bulk states, while the unoccupied edge states (flat band above EF) have lower energies than the unoccupied bulk states.

Fig. 2. The band structures of (a) pristine Z6 × 6, (b) N-Z6 × 6 with N dopant substituted at site 1, and (c) N-Z6 × 6 with N dopant substituted at site 6 in AFM configurations. The black and red lines correspond to spin-up and spin-down components, respectively. Fermi level is denoted by horizontal blue dotted line. (d) Integrated spin density of N-Z6 × 6 with N doped at site 1 in AFM configuration. The spin density is integrated in the normal direction to the ribbon plane, with red and blue colors representing the up and down spin components, respectively. The band structures of (e) pristine Z6 × 6, (f) N-Z6 × 6 with N dopant substituted at site 1, and (g) N-Z6 × 6 with N dopant substituted at site 6 in FM configurations. (h) Integrated spin density of N-Z6 × 6 with N doped at site 6 in FM configuration.

The metallic FM configuration of ZGNRs is less stable than the semiconducting AFM configuration.[4] As shown in Fig. 2(e), two bulk bands cross at EF, whereas two flat bands (edge states) are further away from EF. In other words, in the FM configuration, the occupied edge states (flat band below EF) have lower energies than part of the occupied bulk states, while the unoccupied edge states (flat band above EF) have higher energies than part of the unoccupied bulk states. Due to the occupation of high-energy bulk states, the total energy of the FM configuration is higher than that of the AFM configuration. In this study, we find that the energy sequence of the localized edge states and the delocalized bulk states influences the magnetism of the substitutionally doped ZGNRs, as will be discussed in detail in the following.

To study the magnetic properties of N-doped ZGNRs, we introduce a single N dopant in a Z6 × 6 supercell, corresponding to a doping concentration of 1.39% (or carrier concentration per unit length of 0.68 e/nm). For the pristine Z6 × 6, the AFM configuration is lower in energy than the FM configuration by 26.8 meV per supercell. When N dopant is substituted near the edge of the ribbon, the ground state of N-Z6 × 6 remains the AFM order, the same as the pristine ribbon. The energy of the AFM configuration is found to be lower than the FM configuration by 14.8, 0.5, and 3.5 meV per supercell when N is doped at sites 1, 2, and 3, respectively. Interestingly, when the N dopant is substituted near the center, the FM configuration becomes the ground state, with energy lower than the AFM configuration by 4.6, 3.4, and 5.5 meV per supercell when N is doped at sites 4, 5, and 6, respectively. In other words, the magnetic order of the doped ZGNRs can be controlled by the N substitutional sites, without further modification to the ribbons.

The AFM–FM transition induced by the N substitutional sites is found to be general to ZGNRs. In our study, we consider N substitutional doping in various ZGNRs with different doping concentrations and different ribbon widths: Z6 × N (N = 1, 2, 3, 4, 5, 6, 10, 15) and ZM × 5 (M = 2, 3, 4, 5, 6, 7, 8). If the N doping concentration is not too high (N > 2 and M > 2), all the N-doped ZGNRs experience the AFM–FM transition induced by N substitutional sites. All the doped ZGNR has an AFM ground state when N atom is doped at the edge (site 1). As the substitutional site moves towards the center of the ribbon, the energy difference between the AFM and FM configurations gets smaller, as shown in Fig. 3. When N is doped near the center (sites 4–8), the energy of the FM configuration becomes lower than that of the AFM configuration. The AFM–FM energy difference is about 2–10 meV per dopant when N is substituted near the center of the ribbons.

Fig. 3. The energy difference between the AFM and FM configurations as a function of substitutional sites (a) for N-Z6 × N and B-Z6 × N nanoribbons, and (b) for N-ZM × 5 nanoribbons.

Next, we discuss the mechanism of the AFM–FM transition. Here we take N-Z6 × 6 as an example. N dopant will introduce an extra electron to the ribbon. When N dopant is substituted at the left edge (site 1), the extra electron will mainly occupy the spin-down edge states localized on the left edge, since the spin-up state on the same edge has been fully occupied. Consequently, the originally empty spin-down states localized on the left edge will shift downwards below EF and becomes occupied as shown in Fig. 2(b), leading to a decrease in the band gap compared to pristine Z6 × 6.[34] When N is substituted at the center of the ribbons (site 6), the extra electron will mainly occupy the bulk states. If the doped ribbon remains the AFM coupling, the low-energy edge states have to shift upwards to maintain the empty feature, while the high-energy bulk states have to shift downwards to become occupied, as shown in Fig. 2(c). This band distortion is energetically unfavorable, and the occupation of the high-energy bulk states significantly increases the total energy of the AFM configuration. On the other hand, in the FM configuration, part of the unoccupied bulk states have lower energies than the unoccupied edge states (Fig. 2(e)). The occupation of these low-energy bulk states by the extra electron is energetically more favorable, with EF slightly shifted upwards as shown in Fig. 2(g). As a result, the FM configuration has a lower energy than the AFM configuration when N is substituted at the center. When N is substituted in the intermediate sites, the extra electron would occupy both of the edge and bulk states (not shown). As the doping site moves to the center, the competition between the occupation of edge states and bulk states leads to the AFM–FM transition in the substitutionally doped ZGNRs.

The unstableness of the FM configuration when N is doped near the edge can also be understood by the competition between the edge and bulk states. Suppose that the doped ZGNR has FM order when N is doped near the edge. Since the extra electron will mainly occupy the edge states, the unoccupied high-energy edge states have to shift downwards to become occupied, while the unoccupied low-energy bulk states should shift upwards to remain empty, as shown in Fig. 2(f). This induces a significant distortion in band structure and an increase in total energy. Therefore, the FM configuration is quite unfavorable. For some N-ZM × N ribbons with N ≤ 6, we failed to obtain a metastable FM configuration when N is substituted near the edge (sites 1 and 2), and thus we have some missing points in Fig. 3.

The magnetic orders of the doped ZGNRs can be visualized by the integrated spin density in Figs. 2(d) and 2(h). When N is doped at the left edge, the ground state remains the AFM order, as seen from Fig. 2(d). The total magnetic moment is almost 1.0 μB, mainly contributed by the extra electron localized on the N dopant. The extra electron has an opposite spin and partially reduces the local magnetic moments on the same edge near the dopant.[12,34] For N-Z6 × N ribbons with N ≤ 4 where doping concentration is relatively high, the spin polarization on the doped edge can be almost completely suppressed.[35] However, the N dopant has little effect on the spins of the other edge, due to the localization of the extra electron. Apparently, the edge doping destroys the spin degeneracy of the edge states. The spin anisotropy of the two edge states is predicted to give rise to spin-polarized transmittance channels.[32] On the other hand, when N is doped near the center in FM configuration, the extra electron mainly occupies the bulk states and has little effect on the spins of the two edges, as shown in Fig. 2(h). When N is doped at site 6 in Z6 × 6, the total magnetic moment is 2.3 μB, similar to the value of 2.8 μB for pristine Z6 × 6. Compared to the edge doping, central doping results in larger magnetic moment and possible larger spin currents.

We also notice that there is no AFM–FM transition at very high doping concentrations or in very narrow ribbons. In N-Z6 × 1 and N-Z6 × 2 ribbons with doping concentrations of 8.3% and 4.2% respectively, the ground state is always AFM regardless of N substitutional sites. In these systems, strong N–N interaction would introduce significant band distortion. The edge and bulk states near EF are not well separated in energy. In very narrow ribbons such as Z2 × N, the edge states and bulk states are spacially not well separated. The extra electron from dopant will occupy both the edge and bulk states regardless of the substitutional sites. In these two cases, there is no noticeable competition between the edge and bulk states induced by N doping sites, and thus no AFM–FM transition can be observed. However, too low doping concentration seems unfavorable to stabilize the FM configuration. In Z6 × 10 ribbon, the FM configuration is more stable than the AFM configuration when N is doped at sites 5 and 6, as shown in Fig. 3(a). In a longer Z6 × 15 supercell, the FM configuration is more stable only when N is doped at site 6. This can be ascribed to the large AFM–FM energy difference for pristine ribbons with large N. The energy difference is 68 meV for Z6 × 15 surpercell. One N dopant is difficult to overcome the large energy difference and reverse the relative stability.

The competition between the edge and bulk states can be generalized to hole doping. By introducing one B dopant at different sites in Z6 × N ribbons (N = 5, 6, 10), we observe similar AFM–FM transition depending on the B substitutional sites (see Fig. 3(a)). Unlike N doping, B dopant introduces an extra hole to the ribbon that prefers to occupy high-energy states. When B is doped near the edge of the ribbon, it is energetically favorable for the extra hole to occupy the edge states in the AFM configuration, and thus one edge band shifts upwards above EF (Fig. 4(a)). The local magnetic moments on the nearby edge C atoms are suppressed, as shown in Fig. 4(b). On the other hand, when B is doped near the center, it is energetically favorable for the hole to occupy the bulk states in the FM configuration, and only the bulk states near EF are shifted upwards as shown in Fig. 4(c). The spin density distribution on the edges is similar to the pristine ribbon in FM configuration (Fig. 4(d)). As seen from Fig. 3(a), the AFM–FM energy difference can reach up to 25 meV per B dopant when the dopant is near the center, larger than the N substitution. Therefore, B substitution is more effective to stabilize the FM configuration in ZGNRs.

Fig. 4. (a) The band structure and (b) the integrated spin density of B-Z6 × 6 with B dopant at site 1 in the AFM configuration. (c) The band structure and (d) the integrated spin density of B-Z6 × 6 with B dopant at site 6 in the FM configuration. The gray, green, and white balls represent C, B, and H atoms, respectively.

The AFM–FM transition induced by substitutional doping sites originates from the competition between the edge and bulk states, taking advantage of the properties that the edge states and the bulk states are separated in energy and in space. Previous studies on carrier doping only made use of the property that the edge and bulk states are separated in energy.[25,26] For carrier doping, the injecting carrier will firstly occupy the low-energy edge states in AFM configuration, and then the high-energy bulk states. Only when the carrier doping concentration exceeds a critical point, the occupation on the bulk states in FM configuration becomes more favorable than in the AFM configuration, leading to the AFM–FM transition.[25,26] The critical doping concentration was found to depend on the width of the ribbon.[25,26] This is consistent with our finding that too low doping concentration is difficult to stabilize the FM order. Sawada et al. suggested that carrier doping can be performed by substitutions of B or N atoms. In N = 10 ribbon with doping concentration of 0.41 e/nm (one dopant in a Z10 × 10 supercell), they claimed that the B or N doped ZGNRs have FM ground state, without considering the substitutional sites.[25] Our study suggests that the substitutional sites should be an additional variable to control the magnetic property of ZGNRs. For N-Z10 × 10 ribbon, our calculations show that the FM configuration is indeed the ground state when N is doped at the center. However, the AFM configuration is more favorable than the FM configuration by 22.8 meV when N is doped at the edge.

As seen from Fig. 3, several substitution sites near the center can stabilize the FM configuration. To study the effect of combination of different substitutional sites, we introduce two B dopants in Z6 × 10 ribbon (the same doping concentration as B-Z6 × 5 with one dopant in Z6 × 5 supercell). The two B dopants are substituted at sites 5 and 6, and kept far away from each other. We consider five different geometries and find that the FM configuration is always more stable than the AFM configuration. The energy difference is 6.2–8.9 meV per B dopant, similar to the energy difference of B-Z6 × 5. This indicates that at this doping concentration each B dopant can contribute 5–10 meV to the energy difference as long as the dopants are separated from each other, in spite that they substitute different central sites.

Previous studies showed that the substrate could affect the magnetic ordering of ZGNRs.[36,37] The experimental confirmation of the edge magnetism of pristine ZGNRs was realized much later than its theoretical prediction,[6,7,9] which was ascribed to the unperfect edge orientations and the possible strong edge–substrate interactions.[6] In our previous studies,[38,39] we found that strong edge–substrate interactions with Si substrates can destroy the edge magnetism of both pristine ZGNRs and N-doped ZGNRs, in agreement with other theoretical studies.[36] Au(111) substrate was predicted to interact with ZGNR edges sufficiently weak and preserve the edge magnetism of ZGNRs,[37] which was confirmed by later experiments.[6] Similarly, we speculate a weak substrate interaction should be necessary to maintain the edge magnetism of the doped ZGNRs. Here we take the BN substrate as an example. For pristine Z6 × 6 on BN sheet, the AFM configuration is more stable than the FM configuration by 31.9 meV. When N is doped at the edge, the AFM configuration is still the ground state, with an energy lower than the FM configuration by 21.7 meV. When N is doped at the center, the FM configuration becomes more stable than the AFM configuration by 5.0 meV. Clearly, the weak interaction with BN substrate can preserve the edge magnetism and the relative stability of the AFM and FM configurations in doped ZGNRs.

4. Conclusion

We have systematically studied the electronic and magnetic properties of N- and B-doped ZGNRs, and reveal an AFM–FM transition induced by the substitutional sites. Within appropriate doping concentration, the ground state of the doped ZGNRs is AFM order when the dopant is substituted near the edge of the ribbons, while its ground state becomes FM order when the dopant is substituted near the center. The AFM–FM transition is well explained by the competition between the edge and bulk states. Due to the energy difference of the edge and bulk states, the extra charge from the dopant is energetically more favorable to occupy the low-energy edge states in the AFM configuration when the dopant is doped near the edge. On the other hand, the extra charge is energetically more favorable to occupy the low-energy bulk state in the FM configuration when the dopant is near the center. We also suggest that suitable substrates with weak van der Waals interactions should be necessary to study the magnetism of doped ZGNRs. Our studies provide a guide to synthesize graphene-based nanostructures with stable FM order via a relatively simple doping method, which is highly desirable for the potential applications of graphene-based spintronic devices.

Reference
[1] Makarova T L 2004 Semiconductors 38 615
[2] Wang W L Meng S Kaxiras E 2008 Nano Lett. 8 241
[3] Zhou J Wang Q Sun Q Jena P 2011 Phys. Rev. B 84 081402
[4] Son YW Cohen M L Louie S G 2006 Phys. Rev. Lett. 97 216803
[5] Li X Wang X Zhang L Lee S Dai H 2008 Science 319 1229
[6] Magda G Z Jin X Hagymási I Vancsó P Osváth Z Nemes-Incze P Hwang C Biró L P Tapasztó L 2014 Nature 514 608
[7] Ruffieux P Wang S Yang B Sánchez-Sánchez C Liu J Dienel T Talirz L Shinde P Pignedoli C A Passerone D Dumslaff T Feng X Müllen K Fasel R 2016 Nature 531 489
[8] Nakada K Fujita M Dresselhaus G Dresselhaus M S 1996 Phys. Rev. B 54 17954
[9] Li Y Y Chen M X Weinert M Li L 2014 Nat. Commun. 5 4311
[10] Cervantes-Sodi F Csányi G Piscanec S Ferrari A C 2008 Phys. Rev. B 77 165427
[11] Martins T B da Silva A J R Miwa R H Fazzio A 2008 Nano Lett. 8 2293
[12] Li Y Zhou Z Shen P Chen Z 2009 ACS Nano 3 1952
[13] Biel B Blase X Triozon F Roche S 2009 Phys. Rev. Lett. 102 096803
[14] Zheng X H Wang R N Song L L Dai Z X Wang X L Zeng Z 2009 Appl. Phys. Lett. 95 123109
[15] Zheng X H Wang X L Abtew T A Zeng Z 2010 J. Phys. Chem. C 114 4190
[16] Huang H Gao G Fu H Zheng A Zou F Ding G Yao K 2017 Sci. Rep. 7 3955
[17] Kan EJ Li Z Yang J Hou J G 2008 J. Am. Chem. Soc. 130 4224
[18] Zhang H Meng S Yang H Li L Fu H Ma W Niu C Sun J Gu C 2015 J. Appl. Phys. 117 113902
[19] Nam Y Cho D Lee J Y 2016 J. Phys. Chem. C 120 11237
[20] Wang D Zhang Z Zhu Z Liang B 2014 Sci. Rep. 4 7587
[21] Castro E V Peres N M R Stauber T Silva N A P 2008 Phys. Rev. Lett. 100 186803
[22] Jiang J Turnbull J Lu W Boguslawski P Bernholc J 2012 J. Chem. Phys. 136 014702
[23] Sarmah A Hobza P 2017 RSC Adv. 7 46604
[24] Kan M Zhou J Sun Q Wang Q Kawazoe Y Jena P 2012 Phys. Rev. B 85 155450
[25] Sawada K Ishii F Saito M Okada S Kawai T 2009 Nano Lett. 9 269
[26] Jung J MacDonald A H 2009 Phys. Rev. B 79 235433
[27] Kresse G Hafner J 1993 Phys. Rev. B 47 558
[28] Ceperley D M Alder B J 1980 Phys. Rev. Lett. 45 566
[29] Blöchl P E 1994 Phys. Rev. B 50 17953
[30] Monkhorst H J Pack J D 1976 Phys. Rev. B 13 5188
[31] Pisani L Chan J A Montanari B Harrison N M 2007 Phys. Rev. B 75 064418
[32] Martins T B Miwa R H da Silva A J R Fazzio A 2007 Phys. Rev. Lett. 98 196803
[33] Lieb E H 1989 Phys. Rev. Lett. 62 1201
[34] Cruz-Silva E Barnett Z M Sumpter B G Meunier V 2011 Phys. Rev. B 83 155445
[35] Huang B Liu F Wu J Gu BL Duan W 2008 Phys. Rev. B 77 153411
[36] Zhang Z Guo W 2009 Appl. Phys. Lett. 95 023107
[37] Li Y Zhang W Morgenstern M Mazzarello R 2013 Phys. Rev. Lett. 110 216804
[38] Li J Yang S Li S S 2015 Chin. Phys. Lett. 32 027101
[39] Li J Yang S Li S S 2015 Chin. Phys. Lett. 32 077102